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Fig. 3. An example of weighting functions for N = 3, n 
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Fig. 4. CCF’s form = 1. 


n X N weight elements are the same absolute value 2 / Vn. And the 
others, 2N(N — 1) weight elements with a value of 0, need not be 
implemented in the weighted adders and do not have an effect on 
the amplitude distribution. Therefore, the implementation of the N 
filters is very easy, and the amplitude distribution of every random 
signal xx (mj ) is best approximated to a normal from the viewpoint 
of kurtosis. 

Fig. 3 is an example of the proposed weighting functions for N 
= 3,n = 16. The CCF’s of random signals generated by these 
weighting functions are given in Fig. 4 when m = 1. Since all 
CCF’s are zero atk = 6)(j = 0, +1, +2, °°°), three random 
signals are made uncorrelated by sampling r(j) every 6 shift 
pulses. 


IV. CONCLUSIONS 


A set of N weighting functions for N digital filters has been pro- 
posed based on an even-shift orthogonal sequence. The filter can 
generate N uncorrelated random signals from a single binary ran- 
dom signal. All distribution functions of the generated random sig- 
nals become the same normal distribution. The generation speed is 
a function of N and does not depend on n. Hence, the faster gen- 
eration of the random signals could be obtained. Moreover, the 
implementation of the filters is very easy because all n x N weight 
elements and input signal are binary. 
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A Fast and Accurate Single Frequency Estimator 


STEVEN KAY 


Abstract—A new frequency estimator for a single complex sinusoid 
in complex white Gaussian noise is proposed. The estimator is more 
computationally efficient than the optimal maximum likelihood esti- 
mator yet attains as good performance at moderately high signal-to- 
noise ratios. Also, the estimator is shown to be related to the linear 
prediction estimator. This relationship is exploited to reveal why the 
linear prediction estimator does not attain the Cramer-Rao bound even 
at high signal-to-noise ratios. 


I. INTRODUCTION 


The estimation of the frequency of a single complex sinusoid in 
white Gaussian noise is a problem in signal processing which has 
received much attention. See [1] for a summary. The optimal max- 
imum likelihood estimator (MLE) is well known to be given by the 
location of the peak of a periodogram. This estimator attains the 
Cramer-—Rao lower bound on variance for a high-enough signal-to- 
noise ratio (SNR). In many instances, however, the computation is 
prohibitive even with an FFT implementation, and so simpler 
methods are desirable. In this correspondence we present an ap- 
proach which is strongly motivated by the recent work of Tretter 
[2]. It is shown that the proposed estimator is computationally much 
simpler than the periodogram, yet attains the Cramer-Rao bound 
for moderately high SNR’s. 

In particular, consider the received data to consist of a single 
complex sinusoid in complex white Gaussian noise, or 


0 Pa ee? ae a (1) 


The amplitude A, frequency wo, and phase @ are deterministic but 
unknown constants. It is the frequency wy that we are interested in 
estimating. The amplitude and phase are considered to be nuisance 
parameters. The noise z, is assumed to be a zero mean complex 
white Gaussian process with z, = z,, + jZ2,- Z1;, Z2, are each real 
Gaussian random variables with zero means, variances of a? #2 
(02 is the variance of z,) and uncorrelated with each other. We now 
assume that the SNR, which is A*/o?, is large, allowing the data 
model of (1) to be replaced by an approximate model which will 
form the basis for the proposed estimator. This approximate model 
is [2] 


xX, = Agel eut to + us) (2) 


where u, is zero mean white Gaussian noise with variance 02/2A’. 
Denoting the phase of x, by z2x,, we have finally 


t=0,1,2,°°:,N-1. (3) 


Having obtained (3), Tretter suggested estimating wo and @ using a 
least squares estimator which is equivalent to an MLE. His ap- 
proach provides the insightful result that frequency and phase es- 
timation is equivalent to linear regression of the phase data. The 
only difficulty is that the phase needs to be unwrapped in comput- 
ing these estimates. This unwrapping, besides adding to the com- 
putation, may prove to be difficult at lower SNR’s. In the next 
section, we show how to avoid phase unwrapping but still attain 
the Cramer-Rao bound. Also, the proposed estimator is shown to 
be an improved version of a correlation or linear prediction esti- 
mator previously studied. 


ZX, = Wot +O + uy, 
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II. DERIVATION OF ESTIMATOR 


Assuming that we wish to estimate only the frequency, we can 
avoid phase unwrapping by considering the differenced phase data 


AL = 2%,4) = 2%, (4) 
fort = 0, 1, , N — 2, which becomes, from (3), 
A, = a + U4) — Uy. (5) 


It is clear from (5) that the problem now is to estimate the mean, 
wg, Of a colored Gaussian noise process. The process is actually a 
moving average with coefficients | and —1. The MLE of wo, which 
is equivalent to the minimum variance unbiased estimator for the 
linear model of (5), is found by minimizing [4] 


J=(A — wl)'C"'(A — wl) (6) 


where A= [Ay A, °** Ay_>]’,1 ={1 1.--- 1)’, and Cis 
the (N — 1) X (N — 1) covariance matrix of A,. The solution to 
this problem is well known and is 
lic'A (7) 
gS: 
canes Gl Oa | | 
Also, it can be shown that the variance of this estimator is 


P l 

Var (a) lo: (8) 
It remains only to explicitly evaluate @ and Var (@)). Note that 
@ is unbiased (let A = wol + uw in (7), where [uw], = 4,4, — u, 
fort = 0, 1, 2, , N — 2), and that @» is a Gaussian random 
variable, being a linear function of the data. To evaluate C~', first 
note that A, is a real moving average process with driving noise 
variance a? / 2A’ and coefficients b) = 1, b} = —1. The covariance 
function is thus 


2 2 
oO 0 
c(0) = ep pee) = ¥ 





c(1) = ¢(-1) = os bob) = — 


c(k)=0 {kl =2 


2A 


The covariance matrix takes on the tridiagonal form 





2 il 0 0 res 0 

2 —l 2 -!l 0 cee @ 
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The inverse is well known with the (i, 7 ) element of the (N — 1) 
x (N — 1) matrix being given by [5] 
2A | y 
(Ch, = =| min (i,j) - 4] lsiijjsN-1 (9) 


7 
“~ 


where min (i, j ) denotes the minimum of / and j. After some al- 
gebra, we have that 


N(N? — 1)4 
ro“ = ———— 
602 
N(N? — 1)47 
Ic 'A= ene 2 w,A, 
0. 1=0 
where 
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Hence, we have from (7) that 


Note that £‘,? w, = 1 since &» is an unbiased estimator (at least 
at high SNR). The frequency estimator may further be written by 
using the equivalence 


A, = £241 — 4%, = LXE Kye (10) 
as 
N-2 
Dew, L2XEX | (11) 
with variance which follows from (8) as 
6 
Var (ao) = 7% ’ (12) 
5 N(N? — 1) 
o? 


“ 


But (12) is identical to the Cramer-Rao bound [6]. Additionally, 
the least squares or MLE estimator of Tretter has also been shown 
to attain the Cramer-Rao bound. It is clear then that @ as given 
by (11) and Tretter’s estimator must be identical. In practice, how- 
ever, (11) is to be preferred since it avoids phase unwrapping. To 
verify this equivalence directly, we may rewrite (11) using (10) as 
| N-1 N-| 

ae She 2 2x, (13) 
which is identical to the linear regression estimator of Tretter. That 
(11) and (13) must be the same estimator is also guaranteed by the 
theorem that if an efficient estimator exists (i.e., it attains the Cra- 
mer-Rao bound), then it must be unique [7]. 

The form of the estimator given by (11) is similar to that of a 
previously proposed estimator as will be discussed in Section III. 
It is of interest to note here that w, is a window which is symmetric 
about the point t = N/2 — 1. Some examples of this window are 
given in Fig. 1. As will be discussed shortly, it is this window 
which is responsible for @p attaining the Cramer-Rao bound. If w, 

= 1/N — 1 were chosen, for example, then a soe ue the sam- 
ple mean of the measurements z2x*x,,,, f= 0, 1, sil 
This choice would neglect the colored noise of (5), nics led to 
the need for C' in (6) and ultimately produced w, in (11). In fact, 





for w, = 1/N — 1, the estimator becomes 
1 N-2 
Oy = y=] py LX Xp 
1 N-2 
a ee 
| 
~ Wy i (Z£xy-1 — 2%) 


which although unbiased [see (3)] can be expected to exhibit a large 
variance due to the lack of averaging. It is easily verified that for 
no windowing 


1 


wh 


Var (a) = (14) 


which follows by using (3). The ratio of variances is found from 
(12) and (14) as 


nowindow 


Var (@o) | 


Var( do) | N(N+1) N 


6(N-1) 6 >) 


window 


For large data records, this loss in performance can be substantial. 
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Fig. 1. Frequency estimator window. 


II]. RELATIONSHIP WITH LINEAR PREDICTION ESTIMATES 
The frequency estimates considered in Section II were 





N-2 
Wg = a W, LEX, 4 | (16) 
t= 
where 
3 af — 1 
3N aa: 
= 1 —- | ————— 
mS NPI N 
2 
and 
1 N-2 
o% =F] py) Ady Xpge (17) 


It is possible to find two additional frequency estimators which are 
equivalent to (16) and (17) at high SNR. These are, respectively, 





N-2 
By = 22 WX Xa (18) 
t= 
and 
, N=? 
Wy = 2 Bite (19) 


To form these new estimators, we have interchanged the operations 
of taking the angle and summation. At high SNR, the performance 
of (18) is identical to (16), and that of (19) is identical to (17) as 
shown in the Appendix. It is interesting to note that (19) has been 
proposed by Lank, Reed, and Pollon [3] and later studied by Jack- 
son and Tufts [8] and by Kay [9] as a linear prediction estimator. 
The variance of (19) is given in [3] and is identical to (14). Our 
results show that a windowed linear prediction estimator as given 
by (18) is optimal in that it achieves the Cramer—Rao bound at high 
SNR. At lower SNR, (16) and (18) are different estimators having 
distinctly different performance. The same can be said about (17) 
and (19). Computer simulations which will be described in Section 
IV show, however, that (16) provides the best overall performance. 


IV. COMPUTER SIMULATIONS 


A computer simulation was performed to compare the perfor- 
mance of the four estimators given in (16)-(19). We term (16) as 
the weighted phase averager, (17) as the unweighted phase aver- 
ager, (18) as the weighted linear prediction estimator, and (19) as 
the linear prediction estimator. A data record of N = 24 points was 
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Fig. 2. Performancy of frequency estimators. 


used, and the mean square error versus SNR (both in decibels) de- 
termined. As a basis for comparison, the Cramer—Rao bound, which 
assumes an unbiased estimator, is plotted, and the performance of 
the MLE or location of the peak of a periodogram was also deter- 
mined. The results are shown in Fig. 2 for wo9/2m = 0.05. As 
predicted by theory for high enough SNR, the MLE as well as the 
weighted phase averager and weighted linear prediction estimator 
attain the Cramer-Rao bound. The threshold for the MLE is about 
—1 dB, while that for the weighted phase averager is 6 dB. The 
weighted linear predictor does not appear to exhibit a sharp thresh- 
old but gradually deteriorates in performance with decreasing SNR. 
The unweighted phase averager and linear predictor exhibit the 
same performance at high SNR which as predicted from (15) is 
about 10 log, N/6 = 6 dB below the Cramer-Rao bound. The 
threshold for the unweighted phase averager is about 6 dB, while 
for the linear predictor no threshold is apparent. 


APPENDIX 
INTERCHANGE OF ANGLE AND SUMMATION OPERATORS 


Consider the estimators 


N-2 
OP = Dw, 2x*x,4) 
t=0 
N-2 
OY? = z a WX" X44 


where LN? w, = 1. 
We will show that for high SNR, these two estimators are iden- 
tical. At high SNR using (4) and (5) 


LX X41 = A, = w + Y, 
where v, = u, 41 — U, So that 
Xe) = Me™ 
N-2 N-2 
t=0 t=0 


Now consider the second estimator 
N-2 
& = 2 LY w,(ee) 
t=0 
N-2 
Lem +g 2 wel! 
1= 


N-2 
Wt 2 os wel. 
1=0 
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Assuming | v,| << 1, 


ta 


wy + Zz Qu w,(1 + jv,) 
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N- 
4? 
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N~2 
=wW + Zz (1 rey a vies) 


N ] 


= &) + arctan 2 W,V;, 
1=0 


N-2 
Wo + Qu WU, 
1=0 
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= ah). 
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Performance Analysis of ESPRIT and TAM in 
Determining the Direction of Arrival of Plane 
Waves in Noise 


BHASKAR D. RAO anp K. V. S. HARI 


Abstract—In this correspondence, two subspace based methods, ES- 
PRIT and the Toeplitz Approximation Method (TAM), for estimating 
the direction of arrival (DOA) of plane waves in white noise in the case 
of a linear equispaced sensor array are evaluated: It is shown that the 
least squares version of ESPRIT and TAM result in the same estimate, 
and are statistically equivalent. It is shown that, asymptotically, the 
estimates obtained using Least Squares ESPRIT and Total Least 
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Squares ESPRIT have the same mean squared error. Expressions for 
the asymptotic mean squared error in the estimates of the direction of 
arrival are derived for the methods. Simple closed-form expressions 
are derived for the one and two source case to get further insight. Com- 
puter simulations are provided to substantiate the analysis. 


I. INTRODUCTION 


Subspace based methods for estimating the Direction of Arrival 
(DOA) of plane waves in noise have been developed and studied 
by a number of researchers [1]-[4]. In this correspondence, we 
analyze the performance of two of these methods, ESPRIT [1] and 
the Toeplitz Approximation Method (TAM) [3]. 

MUSIC was the first method that showed the benefits of using a 
subspace based approach [4]. Due to length and technical consid- 
erations, an evaluation of MUSIC and a comparison to the above 
methods is the subject of [5]. Some theoretical results comparing 
MUSIC and the Minimum-Norm method can be found in [6] and 
[7] where a characterization of the methods was done by examining 
the null spectrum. Our work, motivated by these papers, charac- 
terizes the error in the estimate of the direction of arrival directly. 
Some comparisons of MUSIC to ESPRIT can be found in [9]. 
ESPRIT is, like MUSIC, a general approach, and was developed 
to overcome some of the computation and prior information re- 
quirements of MUSIC [1], [8], [9]. Here we only consider its per- 
formance in the context of a linear equispaced sensor array. The 
use of TAM for DOA estimation was first suggested in [3]. A key 
feature of the method is that it is based on a state space model. In 
fact, it will be shown that for the linear equispaced sensor array 
case, ESPRIT can also be described using this formulation. 

The organization of the correspondence is as follows. First, TAM 
and ESPRIT are briefly reviewed and their relationship is estab- 
lished. It is shown that the least squares version of ESPRIT and 
TAM are statistically equivalent. It is then shown that, asymptot- 
ically, the estimates obtained using Least Squares ESPRIT and To- 
tal Least Squares ESPRIT have the same mean squared error. 
Expressions for the asymptotic mean squared error in the estimates 
of the DOA are then derived. The results are specialized for the 
one and two source case leading to interesting insights. Simulation 
results are presented and they support the analysis. 


II. PROBLEM FORMULATION 


The problem of estimating the direction of arrival of M inco- 
herent plane waves incident on a linear equispaced array of L sen- 
sors is considered in this correspondence. For the kth observation 
period (snapshot), the spatial samples of the signal plus noise are 
given by 


T _ (k) (kK) (k) 
tall be eS) al 
M M M 
k ae oe ee 
~ 2 P; 2 a poet es a p} eit oe zy Nis (1) 
i= i= (= 


where w; = 2%(d/X) sin 0;, d being the separation between sen- 
sors, \ the wavelength of the incident signal, and 0; the direction 
of arrival. It is assumed that N, is a mean zero Gaussian random 
vector with independent elements, i.e., E[N,N 7] = o2/. The noise 
is assumed independent of the complex signal amplitudes p ae which 
are also modeled as being jointly Gaussian. The covariance matrix 
P of the amplitudes whose elements are P,, where P; = 
E[ p\p i] is assumed to be of rank M and has distinct eigenval- 
ues. In this correspondence, E[-] and the overbar ‘‘—’’ will be 
used interchangeably to denote the expectation operator. ' 


'Tn this correspondence, T is used to denote transpose, * to denote com- 
plex conjugate, H to denote complex conjugate transpose. Also * is used to 
denote estimates, and the subscript s to denote parameters associated with 
the signal alone. 
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